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A DIGITAL COMPUTER PROGRAM FOR CALCULATING STEADY TEMPERATURE 
AND DENSITY DISTRIBUTIONS IN A GAS CONTAINING HEAT SOURCES 

by Harry W. Davison 
Lewis Research Center 

SUMMARY 

A digital computer program that allows calculation of axial temperature and density 
distributions in a stagnant gas enclosed within an annulus is presented. The program, 
called STADDIG, is written in FORTRAN IV language for the IBM 7094 computer and is 
based on a steady-state one- dimensional (radial) heat-transfer calculation using cylindri- 
cal coordinates. The gas density is related to the gas temperature through the perfect 
gas law. The program allows for conduction across the gas and across the container 
walls. Heat is dissipated from the container walls by forced- convection cooling with an 
incompressible coolant. Heat sources are included in the coolant, the gas, and the con- 
tainer walls. 

The axial distribution of gas temperature and density is determined by successive 
radial heat-transfer calculations at different axial positions. Axial heat conduction is 
neglected. A typical problem containing 20 axial increments requires about 4 seconds of 
running time when the gas densities and thermal conductivities are required to converge 
within 0. 1 and 1 percent, respectively. Changing the convergence requirements has no 
appreciable influence on the problem running time. 


INTRODUCTION 

The determination of temperature and density distributions is a typical problem that 
arises when a neutron-absorbing gas in a nuclear reactor is exposed to a neutron flux. 
The spatial distribution of the mass of gas (gas density) must be known because of its in- 
fluence on the spatial distribution of neutrons and the resulting effects on reactor criti- 
cality. 

This report presents a method of determining the temperature and density distribu- 
tions in a stagnant gas containing distributed heat sources. A digital computer program, 



called STADDIG, allows calculation of temperature and density distributions in a gas ex- 
posed to a neutron flux. The volumetric heating rates in the gas and the container walls 
are expressed as functions of the gas density and the neutron flux. The program 
STADDIG was used in the analysis of a helium 3 control system for a nuclear rocket re- 
actor (ref. 1). 

The mathematical models, including the basic assumptions and a summary of equa- 
tions used in the computer program STADDIG, are discussed. Symbols are defined in 
appendix A; the basic equations are derived in appendix B. The computer program, in- 
cluding details of the input and output data, is discussed in appendix C; and a sample cal- 
culation, which demonstrates the preparation of the input data and shows the format of 
the output, is given in appendix D. Figure 1 contains a flow sheet for the computer pro- 
gram, including references to pertinent tables and equations. 


DESCRIPTION OF MATHEMATICAL MODELS 

The mathematical model provides a tool for calculating axial temperature and density 
distributions in a stagnant gas enclosed within an annulus formed by two concentric tubes, 
as illustrated in figure 2. Heat is generated in the gas and in the tube walls. The heat 
generated in the gas is produced in the exoergic reaction between the neutrons and the 
gas. The energy exists as kinetic energy of the reaction products. (Only two products 
are allowed in STADDIG.) As these reaction products or projectiles undergo collisions 
with the gas, energy is transferred to the gas. All neutron-gas reactions are assumed 
to occur at one neutron energy (that is, to be monoenergetic). The heat generated in the 
tubes has two sources: (1) the absorption of the remainder of the reaction product energy 
and (2) other heat- production reactions, which are assumed to be proportional to the neu- 
tron flux, such as neutron and gamma reactions in the metal. The gas and the tubes are 
cooled by an incompressible fluid flowing inside the inner tube and outside the outer tube. 
One-dimensional (radial) steady- state heat transfer is assumed in calculating a radially 
averaged gas temperature. The axial distribution of the radially averaged gas tempera- 
ture is determined by successive radial heat-transfer calculations at different axial posi- 
tions and includes the heat absorbed by the coolant. The axial gas density distribution, 
for a specified gas pressure, is calculated by the modified perfect-gas relation p = zp&T 
and by the radial gas temperature distribution. 

The model has the following limitations: 

(1) The volumetric heat- generation rates and the thermal conductivities for the gas 
and the metal tubes are independent of radial position. 

(2) The radially averaged gas density is assumed to be inversely proportional to the 
radially averaged gas temperature. 

(3) Heat is transferred through the gas to the tube walls by conduction only but is re- 
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moved from the tube walls by forced convection to an incompressible coolant. 

(4) Heat conduction in the axial direction in the gas and the tubes is neglected. 

The equations outlined below are incorporated in the computer program to calculate 
(1) the radially averaged gas temperature at any axial position, (2) the axial distribution 
of gas density, and (3) the total mass of gas enclosed in an annulus for any selected gas 
pressure. 


Radially Averaged Gas Temperature 

The radial temperature distribution through the gas is determined from the solution 
of the Poisson equation in cylindrical coordinates, with the azimuthal and axial depend- 
ence neglected. The following equations, are required to calculate the radially averaged 
gas temperature. Some of these equations are derived in appendix B. 


T = 



( R 3 - R 2) 
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(See eqs. (BIO) and (B13).) 
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(See eqs. (B21), (B22), and (B23).) 

The thermal conductivity of the gas is approximated by a linear function of the 
average gas temperature T: 


k = a + bT 


( 6 ) 


t k x . = T. , . + 3 j. _ , m > I# j 

b,i,j vApCp 


The equation for 



is presented in appendix B (eq. (B14).) 
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Axial Distribution of Gas Density 

The radially averaged gas density at any position is represented by the relation 
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( 8 ) 


The gas density is related to the gas temperature through the perfect gas law: 

p 

g z£?T 


(9) 


When the value of p^ from equation (9) is substituted into equation (8), the resulting 
integral is difficult to solve in closed form. However, for small radial perturbations in 
gas density or temperature, the average gas density can be approximated by 



z#T 


( 10 ) 
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The average gas temperature T, the average density p , and the heat- generation 
rates q are related through equations (1) to (7), (10), and (B14). These equations are 
iteratively solved for a specific system pressure p to obtain the radially averaged gas 
density. The effect of the axial neutron flux distribution on gas temperature and density 
is transmitted through the heat- generation rates in the gas q^, the tube walls q m , and 
the coolant H. .. The axial distribution of the gas density is obtained by solving the pre- 
ceding system of equations at successive axial positions, with different neutron fluxes but 
the same system pressures. The axial density distribution can be calculated for several 
different system pressures by inserting a pressure increment Apg and the number of 
pressures to be investigated (KOUNT) as input data. Each system pressure Pg is calcu- 
lated from 


pg = pg >s - 4 t>k 


where 


Ap k = ( K “ x ) A Pg 


and 


1 < K < KOUNT 


Total Mass of Gas 

The total mass of gas in a single annulus W is 

W=IY C 


The axially averaged gas density r is 


r = 



(ii) 


( 12 ) 


The annulus volume r is 

c 
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( 13 ) 


T c = ff ( R 3- R D L 
CONCLUDING REMARKS 

A digital computer program called STADDIG written in FORTRAN IV language is 
presented. This program allows calculation of steady- state axial temperature and density 
distributions in a stagnant gas enclosed within an annular space between two concentric 
tubes. Variable heat generation in the axial direction is allowed in both the gas and the 
tubes. Heat is dissipated in the gas by conduction to the tubes. The tubes are cooled by 
forced convection to an incompressible fluid. 

The program includes calculation of coolant temperature distributions, tube wall 
temperatures, friction pressure drop in coolant channels, and total mass of gas in an 
annulus. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, April 25, 1967, 

122-28-02-04-22. 
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APPENDIX A 


SYMBOLS 


o 

A cross-sectional area, cm 

Avogadro number 

a constant used in correlation for 

kg, W/(cm)(°K) 

b constant used in correlation for 

k g , W/(cm)(°K 2 ) 

C arbitrary constant 

C heat capacity, W/(sec)(g)(°K) 

r 

Cg coefficient used in empirical stop- 
ping power - energy correlation 

Cg exponent used in empirical stopping 
power - energy correlation 

D g equivalent diameter, cm 

E energy released after absorbtion of 

neutron, MeV; J 

f friction factor 

G mass flow rate, g/sec 

H heat- generation rate, W 

h forced- convection heat-transfer 

coefficient, W/(cm 2 )(°K) 

K proportionality constant, 

W/(sec)(cm) 

k thermal conductivity, W/(cm)(°K) 

L length of gas annulus, cm 

M molecular weight, g/mole 

Nu Nusselt number 

n number of axial increments 

Pr Prandtl number 


p pressure, N/cm 2 

q volumetric heat- generation rate, 

W/cm 3 

R specific radial field position, cm 

AR gas tube wall thickness, cm 

3 universal gas constant, N(cm)/ 

(g)(mole)(°K) 

Re Reynolds number 

Rjyj radius at maximum temperature, cm 

Rg radius of fuel or, for annuluar cells, 
radius of tube surrounding outer 
tube forming gas annulus, cm 

Rg center -to- center spacing of gas- 

containing tubes, cm 

r general radial field position, cm 

S fraction of E absorbed in gas 

T temperature, °K 

v coolant velocity, cm/sec 

W mass of gas in annulus, g 

WP wetted perimeter, cm 

X axial length, cm 

AX increment of axial length, cm 

z compressibility factor 

/3 linear range of particle, cm 

O 

r axially averaged gas density, g/cm 
jt l viscosity, g/ (cm) (sec) 

O 

p density, g/cm 

0 

p radially averaged density, g/cm 
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a microscopic reaction cross 

2 

section for neutrons, b; cm 

r volume, cm^ 

<p neutron flux, neutrons/ 

o 

(cm )(sec) 

2 

if/ range of particle, g/cm 

Subscripts: 

a portion of heating rate in tube 

walls due to absorption of 
reaction products 

b bulk water conditions 

c calculated values of density 

and conductivity 

f average conditions in convec- 

tive film 

g average conditions in gas 


I inner coolant channel 

i estimated 

j axial location 

k pressure increment 

M maximum 

m tube material 

o outer coolant channel 

p, t reaction products produced in 

gas-neutron reaction 

s initial condition 

1, 2, 3, 4 inner and outer surfaces of 

inner and outer tubes form- 
ing annulus 

Y portion of heating rate in tube 

walls that is proportional 
to cp 
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APPENDIX B 


DERIVATION OF EQUATIONS 

The equations used to calculate the radially averaged gas temperature in STADDIG 
are derived herein. The one-dimensional Poisson equation is used to calculate the ra- 
dial temperature distribution through a stagnant gas containing heat sources and held in 
an annulus formed by two concentric tubes (fig. 2). 


+ 1 dT + q = 0 R < r < R* (Bl) 

dr 2 rdr k 


T = -^ — + C. In r + C 9 
k 4 1 * 


(B2) 


The boundary conditions are 


T = T 2 at r = R 2 

— = 0 at r = R m 

dr 


where R^ is the radius at maximum temperature. Inherent in the second boundary con- 
dition is the assumption that the radius of maximum temperature actually occurs in the 
gas (R 3 > R m > R 2 ). 

For certain coolant-flow distributions, the maximum temperature could occur in the 
wall of either tube rather than in the gas. Although the equations governing this situation 
are not derived herein, they are included in the computer program. The equations can 
be easily derived by the techniques presented herein. 

The gas temperature at any radial position is 


T 


= t 2 + 


*g 


4k 


g 



The average gas temperature in the annulus is 


(B3) 
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T 



T 






(B4) 


(B5) 


2 

The variables R^, T 2 , qg, and must be determined. All other variables in 
equation (B5) are given in the statement of the problem to be solved. 

The maximum temperature T M is calculated as a function of R^j and T^ j and as 
a function of R^j and T^ Q . The two relations for Tj^ are equated, and the resulting 

2 , 

equation is solved for R^j. 


T M - T b,I + - T 2> + ^2 - T l> + <*1 - T b,I> 


' T b,o + < T M- T 3> + < T 3- T 4> + ( T 4- T 


b.o'J 


(B6) 


The term T M - T 2 is determined by equation (B3). 


(Tm - T 2> - 


= _li. 

4k 

S 


' r 2-4 +2r L ln ^l 

ti 2 l 


(B7) 


The term T^ - Tg is determined by using the following boundary conditions with 
equation (B2): 


T = Tg at r = Rg 


dT 

dr 


= 0 


at r = R 


M 
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(T M ■ T 3 ) = ^ ( R 3 ' R M + 2R M ln “] 

4k g \ R 3y 


(B8) 


The following boundary conditions were used with equation (B2) to solve for Tg - T^: 


T = T, at R = R, 


■ 

m dr 2R 0 


at r = R 0 


The second boundary condition represents constant heat flux across the boundary be- 
tween the gas and the inner tube. 


T ' T i + i^( R ?- r2 ) 


«- R 2> S + ^m K 


2k 
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( R h - glK + 9 m R 
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Rrt 
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(BIO) 


The temperature difference Tg - T 4 is determined from 


T = T 4 at r = R 4 


k ™ — 
m dr 


dT _ ( R M - R 3) q ! 


g 


2R, 


at r = R 


> (R 3 < r < R 4 ) 


The second boundary condition represents constant heat flux across the boundary between 
the gas and the outer tube. 



(Bll) 
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( t 3 - T 4 
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- R l)<lg «m : R t 


2k m 
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R 3 

In — 
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(B12) 


The temperature rise across the outer and inner convective coolant films T 4 ~ 0 


and T, - T, T are determined by a heat balance: 
1 b, I 


< T 4' T b,o> = 


_ ( R l - *&)% 4 H - R IK 


2R 4 h b,o 


(Tx- T b,I> = 


. « - «ik * ( R l - R !)% 


2R l h b, I 


(B13) 


Substituting equations (B7), (B8), (BIO), (B12), and (B13) into equation (B6) and solv- 
ing for yield 


R 2 - 
r M " 


q g 
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/ R 3 \ k e , 

lnl — ] + — s. In 
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where 
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Forced-Convection Film Coefficient 


The Colburn equation (ref. 2) was used to approximate the forced- convection film 
coefficients for both the inside and outside coolant channels: 

Nu f = 0.023 Pr 1//3 Re 0 * 8 

h = 0- 023 Pr 1//3 Re°- 8 k ^ B15 ^ 

D e 

Volumetric Heat-Generation Rate in Gas 

When the gas absorbs a neutron, an unstable intermediate product may be formed. 
This unstable intermediate product decomposes into two secondary recoil products (only 
two recoil products are allowed in the computer program). Heat is generated in the gas 
as the result of the slowing down of the recoil products. The volumetric heat- generation 
rate is determined from (1) the calculated reaction rate between the gas and the neu- 
trons, (2) the energy released in each reaction, and (3) the fraction of energy absorbed 
in the gas. 


q 


g 


-13 Ej ^N CTP £® t/7 
1.602x10 16 5L_I — 

M 


(B 16) 


The values of E, a, cp, and M must be specified as input data for the calculation. The 
fraction of the energy absorbed in the gas S may be either specified as input data or cal- 
culated. 

If S is to be calculated, the program incorporates the model proposed in refer- 
ence 1. In this model it is assumed that the average distance a projectile travels j3 in 
an annulus is the radius of a sphere whose volume is equal to the volume in which a pro- 
jectile can be stopped in an annulus: 


P - 



P 


g 


(B17) 


The energy lost by the reaction products while slowing down in the gas is determined 
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by integrating an empirical approximation for the stopping powers dE/d(Pg/3) of each 
projectile. The following stopping- power approximation is used in the program: 

^ = -CjjPgE 06 (B18) 

d/3 

Integration of equation (B18) yields 

1-C R 1-C R 

E s - E b = C 5 (1 - C 6 )p g /3 (B19) 

where E„ is the initial energy of the projectile, 
s 

The fraction S of the energy deposited in the gas is 


S = 


< E s - E >p + - E >t 

E s,P + V 


(B20) 


The subscripts p and t refer to the two reaction products. The initial energy of the 
projectiles is specified as input data. The final energies of the projectiles are obtained 
from equation (B19). 


Volumetric Heat-Generation Rate in Tube Walls 


Heat is generated in the tube walls and is dissipated at the wall surface by forced 
convection. Two types of volumetric heat sources are provided. The first q^ is pro- 
portional to neutron flux (e. g. , heating due to absorbtion of primary gamma photons): 

<iy - K<? (B21) 

The second q results from absorption of the portion of the reaction products or pro- 

a A 

jectiles that is not absorbed in the gas: 


«a 


yi - siK - r|) 

K r ! - «i + r I - 


(B22) 
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where is obtained from equation (B16) and S is obtained from equation (B20). 

The total volumetric heating rate in the tube walls is 

q m = + q a (B23) 

Coolant Temperature 

The axial distribution of coolant temperature must be determined before the gas tem- 
perature and density can be determined. The coolant temperature at each elevation is cal- 
culated from the preceding coolant temperature, the power generation at each section, 
and the coolant flow rate: 


T _ T , H b,j +H g,j +H m,i 

T b,j- T b,(j-1) + — 


(B24) 


GC 


p,b 


The heat generated in the coolant is 


H b, j = q b A o AX I7 when r > R 4 


<P 


H b j = q b A I AX when r < 


<P 


(B25) 


where <p/cp is the relative neutron flux at any axial location. 

At any axial location j, the heat- generation rates in the gas in the regions inside 
(H . ,) and outside (H . ) the radius of maximum temperature are 




and 


H g,j,I = ’' q g iX ( R M' R l) whetl R 2 <r<R M 


H g, i, o = ’"‘g aX ( R 3 - “m) when H 3 > r > 


(B26) 


The heat- generation rates in the inner and outer tubes are 
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and 


H m,j - "1m iX ( R 2 ' R l) when R 1 < r < R 2 


V j ' " q m AX ( R 4 - R i)) when R 4 > r > R 3 


(B27) 


The coolant physical properties are assumed to be constant along the length of the tubes. 


Equivalent Diameters of Coolant Channels 

The equivalent diameter D g of a channel is defined as four times the cross- 
sectional flow area A, divided by the wetted perimeter WP. The equivalent diameter 
D g of the inner channel is the diameter of the inner tube. The equivalent diameter of 
the outer channel depends on the lattice geometry in which the tube is located. Any of 
three separate geometries is permitted in this program: square, equilateral triangular, 
or annular. These geometries are illustrated in figure 3. For the first two geometries, 
the gas annulus is located at the center of the cell surrounded by other cylindrical tubes, 
such as fuel elements. In the annular geometry, the outer tube of the gas annulus is con- 
centric with the larger tube that surrounds it. Equations for cross-sectional flow area 
A, wetted perimeter WP, and equivalent diameter D e are given in table I. 

The lattice geometry is selected in the computer program by inserting the appropri- 
ate value of RG, as indicated in appendix C. 


TABLE I. - TUBE GEOMETRY 


Geometry 


Equation for 

- 

Figure 


Cross-sectional area, 

A 2 

cm 

Wetted 

perimeter, 

WP, 

cm 

Equivalent diameter, 
cm 


Equilateral 

triangle 

1.299R2-I (r2 +2R 2) 

+ 2R 4 ) 

5. 196 Rg - 2jt(r| + 2R^) 

4(a) 

77 (R g + 2R 4 ) 

Square 

R 6 ' *( R 5 + R4) 

2?7(R 5 + R 4 ) 

2R 6- 2*(r2 + r 2) 

4(b) 

Jt(R 5 + R 4 ) 

Annular 

77 ( R 5 ' **4) 

2fir(Rg + R 4 ) 

2(Rg - R 4 ) 

4(c) 
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Friction Pressure Drop in Coolant Channels 


The irreversible pressure loss in the inner and outer coolant channels was calculated 
by the Fanning equation for isothermal flow (ref. 3) : 


Ap = 


2fLy2p b 


(B28) 


v 


G 

P b A 


(B29) 


The friction factor 


f is approximated by 


f = 0. 046 Re -0- 2 (B30) 

where Re is the Reynolds number. This friction- factor correlation is valid for fluids 
flowing in smooth tubes and for Reynolds numbers between 10 4 and 10®. 
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APPENDIX C 


DESCRIPTION OF COMPUTER PROGRAM 


The axial distribution of gas temperatures and densities and the total mass of gas in 
an annulus are calculated by equations (10) to (13) and (B14). The only iterative proce- 
dures in the computer program are associated with the determination of the radially aver- 
aged gas density and the thermal conductivity of the gas. The initial estimate of the gas 
density is 1 milligram per cubic centimeter. All density dependent variables, including 
average gas temperature, are calculated by using this initial guess. The density p is 
calculated (eq. (10)) from the temperature, the pressure, the molecular weight, and the 
compressibility of the gas and is compared with the estimated value. The absolute values 
of the fractional difference between estimated and calculated values of density are com- 
pared with a convergence parameter CONVl which is part of the input data. When 



> CONVl 


a new value of the density is estimated and the calculations are repeated. If the result 
is equal to or less than the convergence parameter, the solution is accepted, and a cal- 
culation is initiated for the next axial increment. 

The initial estimate of the thermal conductivity of the gas k . is supplied as input 
data. This value is used to calculate the average gas temperature at any axial location. 

A new thermal conductivity k is calculated by using equation (6) and the average gas 
temperature. The absolute value of the fractional difference between calculated and esti- 
mated conductivity is compared with another convergence parameter CONV2, which is 
part of the input data. When 



> CONV2 


a new value of k . is estimated, and the calculation is repeated until 

< CONV2 


k g, c ~ k g, i 


k g,i 


- 1 - 5 

The values of CONVl and CONV2 were varied between 10 and 10 to determine 
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their effect on problem running time. The maximum variation in problem running time, 
when only the convergence parameters were varied, was 3 to 5 seconds. 

The input data for each problem require no fewer than 10 or 12 cards, depending 
on the option chosen for calculating the fraction of the total reaction {Product energy de- 
posited in the gas. The input and output variables and options and the input data format 
are described in the two lists that follow. The character immediately following the 
FORTRAN variable name is the algebraic character used in the preceeding equations. 


Input 


Card 

1 

FORTRAN 

variable 

Algebraic 

character 

Format 

72H 

Description 
Problem title 

Units 

2 

WGAM 

% 

F10. 5 

Average heating rate 
in coolant 

W/cm 3 


FLOW1 

G I 

F10. 5 

Coolant mass flow 
rate to inner channel 

g/sec 


FLOW2 

G o 

F10. 5 

Coolant mass flow 
rate to outer channel 

g/sec 


DIST 

L 

F10. 5 

Total length of gas 
annulus 

cm 


CONM 

k 

m 

F10. 5 

Thermal conductivity 
of tube walls 

W/(cm)(°K) 


CONG 

k 

S 

F10. 5 

Thermal conductivity 
of gas 

W/ (cm) (°K) 

3 

RAD2 

r 2 

F10. 5 

Inner radius of an- 
nulus 

cm 


RAD3 

R 3 

F10. 5 

Outer radius of an- 
nulus 

cm 


DELR 

AR 

F10. 5 

Thickness of tube 
walls 

cm 


PRESS 

P g 

Ap 

F10. 5 

Gas system pressure 

N/cm 3 


DELP 

F10. 5 

Pressure increment 
to be investigated 

N/cm 3 
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Card 


4 


FORTRAN 

variable 

Algebraic 

character 

Format 

Description 

Units 

RG 


F10. 5 

Option chosen to se- 
lect proper lattice 
geometry 

Options 

RG Geometry 

1. 0 square 

2. 0 annular 

3. 0 triangular 


SIGMAG 

a 

F10. 5 

Microscopic reaction 
cross section for 
neutrons 

b 

ATOMWT 

M g 

F10. 5 

Molecular weight of 
gas 

g/mole 

TB1 

T b,i 

F10. 5 

Coolant inlet temper- 
ature for inner chan- 
nel 

°K 

TB2 

t. 

D, O 

F10. 5 

Coolant inlet temper- 
ature for outer chan- 
nel 

°K 

CODE 

— 

F10. 5 

Constant which must 
be 1.0 


EPSOG 

E 

F10. 5 

Kinetic energy of re- 
action recoil parti- 
cles 

MeV/reaction 

RAD 5 

R 5 

F10. 5 

Outer radius of cy- 
lindrical element 
(fig. 3) or, for an- 
nular cell geometry, 
inner radius of tube 
surrounding gas- 
filled tube 

cm 

RAD6 

R 6 

F10. 5 

Center-to- center 
spacing of gas- 
containing tubes 

cm 
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Card 


FORTRAN 

variable 


Algebraic Format 
character 


Description 


Units 


6 


7 


SPLIT S F10. 5 


CHOICE F10. 5 


CON5 q m E10. 5 

AFLUX 7p E10. 5 

IPEAK 13 

MM n 13 

KOUNT 13 


Fraction of total en- 

ergy produced in 
neutron- gas reaction 
that is absorbed in 
gas 

Option chosen to se- 

lect method of deter- 
mining SPLIT 

Options 

CHOICE Method of 
determin- 
ing SPLIT 

1. 0 use value 

supplied in 
card 5 
(SPLIT) 

0. 0 calculate 

using 
model de- 
scribed 
previously 

Average volumetric W/cm^ 

heat- generation rate 
in tube walls 

Average thermal neu- neutrons . 

O 

tr on flux (cm ) (sec) 

Estimated position of 

maximum flux (any 
number between 1 
and MM) 

Number of axial in- 

crements to be inves- 
tigated 

Number of pressure 

values to be inves- 
tigated 
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Card 

FORTRAN 

variable 

Algebraic 

character 

Format 

Description 

Units 

8 

RHOC 

p b 

F10. 5 

Coolant density 

g/cm 3 


vise 

Mv, 

F10. 5 

Coolant viscosity 

g/(cm)(sec) 


CONDC 

U 

k b 

F10. 5 

Coolant thermal con- 
ductivity 

W/(cm)(°K) 


CAPC 


F10. 5 

Coolant heat capacity 

(W)(sec)/(g)(°K) 

9 

ACON 

F> u 
a 

E10. 5 

Empirical constant 
for linear fit of gas 
thermal conductivity 

W/(cm)(°K) 


BCON 

b 

E10. 5 

Empirical coefficient 
for linear fit of gas 
thermal conductivity 

W/(cm) ^K 2 ) 


ZCON 

z 

E10. 5 

Compressibility fac- 
tor for gas 



CONV1 

— 

E10. 5 

Convergence param- 
eter for gas density 



CONV2 


E10. 5 

Convergence param- 
eter for gas thermal 
conductivity 


a 10 

BFLUX (I) 

<P 

(2 OX, 
E16.8) 

Relative thermal 
neutron flux at axial 



position j 


The following two cards must be supplied when SPLIT is to be calculated (CHOICE 
0. 0), but must be omitted when CHOICE = 1.0. 
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RAEXPP 


RAEXPT 


'6, p 


'6,t 


E10. 4 Exponent used in em- 
pirical correlation 
between stopping 
power and range for 
particle p 

E10. 4 Exponent used in em- 
pirical correlation 
between stopping 
power and range for 
particle t 


a One card is required for each value of thermal neutron flux. 
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I 


Values for MM of the flux must be provided. 


Card FORTRAN Algebraic Format Description Units 



variable 

character 





RACOFP 

C 5,P 

E10. 4 

Coefficient used in 
empirical correla- 
tion between stopping 
power and range for 
particle p 



RACOFT 

C 5,t 

E10.4 

Coefficient used in 
empirical correla- 
tion between stopping 
power and range for 
particle t 

9 

12 

RANGP 


E10.4 

Range of particle p 

g/cm Z 

9 


RANGT 

*t 

E10. 4 

Range of particle t 

g/crn 


EOPRO 

E s, p 

E10.4 

Initial kinetic energy 
of particle p 

MeV 


EOTRI 

e m 

E10.4 

Output 

Initial kinetic energy 
of particle t 

MeV 

Line 

FORTRAN 

variable 

Algebraic 

character 

Format 

Description 

Units 


1 to 3 

List of input data from cards 1 to 3 

4 

List of input data from card 4 

Length of each axial increment in centimeters (DELX) 

5 

List of input data from card 5 

6 CON4 Average energy pro 

duction per unit of 
volume per unit of 
gas density due to 
neutron interaction 
with gas 


(W/cm 1 * 3 * 5 6 ) 

(g/cm 3 ) 
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II II 


II III 


I III 


Card FORTRAN Algebraic Format Description 

variable character 

List of input data from card 6 
7 to 9 

List of input data from cards 7 to 9 
b 10 to 12 

List of input fluxes from cards 10 to MM 
13 and 14 

List of input data from cards 13 and 14 
b 15 to 17 


List of input fluxes normalized to average neutron flux 


18 

MASS 

W 

Total mass of gas in 

annulus 

19 

PRESS 

P g 

Total gas pressure 


AVG DENSITY 

r 

Average gas density 


in annulus 

Rg Inner radius of an- 

nulus 

Ratio of maximum 

gas density in an- 
nulus to minimum 
gas density 

Vj Coolant velocity in 

inner coolant channel 

V 0 Coolant velocity in 

outer coolant channel 

hj Forced- convective 

film coefficient in 
inner channel 

h Q Forced- convective 

film coefficient in 
outer channel 

Total volumetric 

heat- generation rate 
of metal in stage MM 

b Greater or fewer lines may be required depending on number of axial 
(MM) required. 


RAD2 

GRADIENT 

20 VEL(IN) 
VEL(OUT) 
FILMH1 

FILMH2 

21 Q2 (METAL) 


Units 


g 

N/cm 2 

g/cm 3 

cm 


cm/sec 
cm/sec 
W/(cm 2 ) (°K) 

W/ (cm 2 ) (°K) 

W/cm 3 

increment 
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Card 


23 


Tabular 


FORTRAN Algebraic 

variable character 

QGAM 

QP+T 

Q1 (INGAS) 

IMINRO 

IMAXRO 

PDRP1 Apj 

PDRP2 Ap 0 

REYIN Rej 

REYOT Re Q 


Format Description 

Portion of 

Q2 (METAL) not due 
to attenuation of re- 
coil particles in 
stage MM 

Portion of 

Q2(METAL) due to 
attentuation of recoil 
particles in stage 
MM 

Volumetric heat- 

generation rate of 
gas in stage MM 

Axial location of 

minimum gas density 

Axial location of 

maximum gas den- 
sity 

Irreversible (fric- 
tion) pressure loss 
in inner coolant 
channel 

Irreversible (fric- 
tion) pressure loss 
in outer channel 

Reynolds number in 

inner channel 

Reynolds number in 

outer channel 


listing: 

ELEVATION Axial distance 

DENSITY p Average gas density 

5 

AVG TEMP Tg Average gas temper- 

ature 


Units 

W/cm 3 

W/cra 3 

W/cm 3 


N/cm 3 

N/cm 3 


cm 

g/cm 3 

°K 



FORTRAN 

variable 

Algebraic 

character 

Format Description 

Units 

MAX TEMP 

T M 

Maximum gas tern- 

perature 

°K 

RMAX 

r m 

Radius at maximum 

temperature 

cm 

TBULK(IN) 

T b,I 

Bulk coolant temper- 

ature in inner cool- 
ant channel 

°K 

TBULK(OUT) 

T b, o 

Bulk coolant temper- 

ature in outer cool- 
ant channel 

°K 

TWALLM1 


Maximum tempera- 

ture of inner tube 
wall 

°K 

TWALLM2 


Maximum tempera- 

ture of outer tube 
wall 

°K 
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APPENDIX D 


I 


SAMPLE CALCULATION 

A sample problem and its solution are presented to demonstrate the preparation of 
input data for use in STADDIG and to illustrate the program output format. Because, in 
this sample problem, preparation of data and output format are emphasized, the accuracy 
of the data and the assumptions is of secondary importance. Therefore, the results 
should not be construed as an exact solution of the problem presented. 


Problem 

A nuclear reactor is controlled with helium 3 as a neutron poison. The helium 3 is 
contained in an annulus between two concentric aluminum tubes (fig. 4) which are cooled 
by water flowing inside the inner tube and outside the outer tube. Heat is generated from 
the 2He^(n, p)^T^ reaction at a rate of 0. 765 MeV (1. 22x10“^ J) per reaction. The 
sample problem is to determine the axial temperature and density distribution of helium 3 
and the total mass of helium 3 held in an annulus, when the system is pressurized to 
41 newtons per square centimeter. 

15 

The average thermal neutron flux in the reactor is 1. 5x10 neutrons per square 
centimeter per second, and the axial flux distribution is shown in figure 5. The neutrons 
are assumed to be monoenergetic; therefore, all reactions are thermal. The tube dimen- 
sions and the properties of helium 3, aluminum, and water are shown in table n. The 
average heating rates in the aluminum and water are 121 and 150 watts per cubic centi- 
meter, respectively. The water mass flow rates to the inner and outer channels are 
190 and 400 grams per second, respectively, and the inlet water temperature is 365° K. 
The range-energy relations for the protons and tritons were approximated by using the 
empirical relation in the program. The values of Cg and Cg were obtained empirically 
by plotting the range- energy data for protons in helium 4 (ref. 4) and by correcting the 
range using the Bragg- Kleeman rule (ref. 5). The same constants were used for the pro- 
tons and the tritons. The values used for the coefficients and the exponents in the range- 
energy correlation were 276 and -0.75, respectively. 

The cylindrical fuel assemblies have an outer diameter of 7. 0 centimeters and are 
placed in a triangular pattern. The gas- containing tubes are located in the interstices 
between fuel elements, as shown in figure 4. The center-to-center spacing of the tubes 
is 4. 6 centimeters. 
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TABLE n. - DATA FOR SAMPLE PROBLEM 


Tube dimensions, cm 

Outer -tube inside diameter 1. 940 

Inner -tube outside diameter 1. 740 

Wall thickness 0.15 

Tube length ‘ 107 

Helium 3 properties 

Microscopic absorbtion cross section, a b 5400 

Molecular weight, g 3.0 

Thermal conductivity, ^ a + bT 

a, W/(cm)(°K) 0.865X10" 3 

b, W/(cm)(°K 2 ) 0. 321X10" 5 

Compressibility factor 0 1.0 

Range of proton, d g/cm 2 0. 755x10" 3 

Range of triton, ^ g/cm 2 0. 866x10"^ 

Initial kinetic energy of proton, e MeV (J) 0. 573 (0. 914x10" ^ 3 ) 

Initial kinetic energy of triton, e MeV (J) 0. 191 (0. 306x10 * 3 ) 

Aluminum properties 

Thermal conductivity, * W/(cm)(°K) 1. 989 

Water proper it es* 

Density, g/cm 3 0. 948 

Viscosity, g/(cm)(sec) 0.26xl0" 2 

Thermal conductivity, W/(cm)(°K) 0.68xl0" 2 

Heat capacity, (W)(sec)/(g)(°K) 4. 186 


a Ref. 4. 

b Values of constants are based on linear fit of corrected conductivity data 
from ref. 2. Conductivity of helium isotopes was assumed to be proportional 
to square root of their atomic weights, as proposed in ref. 5. 
c Ref. 6. 
d Ref. 7. 
e Ref. 1. 
f Ref. 2. 


Solution 

The variables given in the statement of the problem are presented in the appropriate 
system of units for use in the computer program. The input data are shown in table EL 
The options, RG = 3. 0 and CHOICE = 0. 0, were selected to permit (1) selection of trian- 
gular spacing and (2) calculation of the fraction of energy deposited in the gas. The con- 
vergence criterion selected for the radially averaged gas density and the thermal conduc- 
tivity of the gas are 0. 1 and 1 percent, respectively. 

The computer output for the problem is given in table IV. The average and maximum 
gas temperatures, the inner and outer tube temperatures, and the coolant temperatures 


28 



TABLE HI. - INPUT FORMAT FOR IBM 7094 COMPUTER PROGRAM 


PROJECT NUMBER ANALYST 


SAMPLE PROBLEM INPUT DATA 


FORTRAN STATEMENT 


STEADY 


15 0.0 

(_-) 1 v--. 1 

0.870 

— I— H 1 1 i 

54 00 . 0 
0 .7 6 4 

— i — i — p— i 1 — 

0,. , 1,3,1 2 E.+ , 

10 21 

1 1 j 

0. 948 3 

— I \ 1 1 1 

0. 8 6 5 0 E- 


T.E.MP ER.A.T 


_ 190, -0 

0 . 97 0 [ 
5. 0 
s'. 5 0 ' 


A.TURE. AND. DENSITY DISTRIBUTION. I. 


40 0.0 
0^ 15 0 | | 
56 5. 0 | 

V. 6 0 


10 7.0 
41. 0 | 
565 .0 
o' - 3 7 ' 


' 1. 9 6 9 
1.5 t 
1^ 0 [ 

V. cT ! 


0.0026 

3. 0 


0 , 5 , , , 1 ,. , 5 , 5 ,+ 1 , 5 , , , , , , , , 

1 . 
0 ■ 0 , 0 , 2 6 , 0 . 0 0 6 8 | 

03 0.321E-05 1 . OOo'eV 

— ! 1 1 1 — l ! 1 1 1 1 — I — ( 1 1 — t — I 1 

0 . 


. 4.1860 

I- 00 l i .000E-05 | 1. PO PE- 02 
,7869E+0p 
. 725 8E+00 

-I 1 1 1 ) V 1 1 1 1 1 1 — I i , i- \ , . 

. 8 5 28E+00 | 

. 9 5 45 E + 0 0 

— * 1 1 1 1 1 I 1 \ 1 1 * ■ i * 1 1 1 1 H 

, 05 28 E+0 0 

—I 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 (- 

. 13 3 7 E+0 0 

—l — t i 1- 4- l — ( 1 — I — l i ! 1 1 — i — i — i 1 H 

, 2017E+00 

_| .+H 1 •—*- | 1 1 1 1 1 1 1 1 1 1 1 f- 

, 25 1 3 E + 0 0 

I * 1 1 i i > * » — ■■+ - < i * * 1 - »- ■ + ■ 

.28 4 6 E+0 0 
. 29 9 3 E+ 0 0 

— « 1 1 »■ t — 4 . - + ■-■- 1 1 , 1 j —4 » 1 1 1 ( f- 

.29 65 E + 0 0 

— I 1 1 1 1 1 1 1 1 ! 4 1 4—1 1 1 1 1 1- 

. 2753E + 00 


. 25 6 OE+ OO 
. 17 7 9E + 0 0 
. 10l3 7 E+0 0 


6 7 8 9 10 li 12 

13 14 15 16 17 18 

19 20 21 22 23 24 

25 26 27 20 29 36 

31 32 33 34 35 36 

37 38 39 4 0 41 42 

A3 UU A5 U6 47 A8 

A9 50 Si 52 53 5A 

55 56 57 58 59 60 

61 62 6364 65 66 

67 68 69 70 71 72 

73 74 75 76 77 78 79 80 









CO 

o 


TABLE m. - Concluded. INPUT FORMAT FOR IBM 7094 COMPUTER PROGRAM 







TABLE IV. - INPUT AND OUTPUT FOR SAMPLE PROBLEM 


INPUT RATA 

STEADY TFMPFRATIJRF AND DENSITY DISTRIBUTION IN A GAS WITH HEAT SOURCES 

WGAMM A* 150*00 FL()W1= 1<50.0 FL0W2* 400.0 DIST* 107.0 CDNM* 1.989 CONG* 0.0026 

RAD?* 0.8700 RA03* 0.9700 DEL R= 0.150C00 PRESS* 41.0 DELP= 1.5 RG= 3.0 

STGMAG* 5400. CG A TDMWT = 3.00 TBi* 365.0 T R? = 365.0 C0DE=1.0 DELX* 5.10 

EPSDG * 0.7640 RAD5* 3.5000 RAD6* 4.6000 SPLIT* 0.3700 CHOICE* 0. 

C0N4=0. l q 91 F 06 CDN5=0.1212F 03 AF LUX* . 15000 E 16 

I PEAK* 10 MU OF INCREMENTS* 21 KOUNT* 1 

RHOC = 0.9483 VISC* 0.0026 CONOC= 0.0068 CAPC* 4.1860 

AC0N=0 .R6 5 OF- 03 BCON=0. 32 1 0F-05 ZCGi\i=0. 1000 E 01 CDNV 1*0 . 10 O0E-O2 CQNV2=0. 1000E-01 

FLUX* 0.7869 0.7758 0.8578 0.9545 1.0528 1.1337 1.2017 1.2513 1.2846 1.2993 

1.7965 1.2753 1.7360 1.1779 1.1037 1.0130 0.9090 0.7890 0.6629 0.5093 

0.4858 

RAFXPP* -0.7500F 00 RAEXPT* -0.7500E 00 RACOFP* 0.27o0£ 03 RACOFT* 0.2760E 03 

RANGP=0*7550t-03 R ANG T*0. 8660F-04 E0PR0=0 . 5730F 00 EOTR I =0 . 1910E 00 

***CALCULATIONS*** 

NORMAL IZFD FLUX 

FLUX* 0.7868 0.7757 0.8527 0.9545 1.0527 1.1336 1.2016 1.2512 1.2845 1.2992 

1.7963 1.2751 1.7359 1.1778 1.1036 1.0129 0.9089 0.7890 0.6628 0.5093 

0.4858 

MASS* 0.071673 

PRESS* 41.00 A VG.DENSITY*^50407E-03 RAD2* J. 870000 GRAD I ENT* 1 . 1 1 46 1 

VEL(IN)* 123.0 VF1 (OUT)* 9d.O FILMH1* 0.90678 FILMH2* 0.82076 

07( METAL )=679521F 07 0GAM*588787F 02 OP* T=907335F 01 Ql UNGAS )=616990E 01 

IMINRO* 14 IMAXRO* 1 PQRP1*.107E 00 PDRP2=.U6E 00 

RFY IN* • 65649F 05 RE YO T=. 346 73E 05 


00 
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TABLE IV - Concluded. INPUT AND OUTPUT FOR SAMPLE PROBLEM 


AXIAL DISTRIBUTION OF RMAX .TEMPERATURES AND DENSITIES 


ELFVATTON 

CMS 

DFN S IT> 
GMS/CC 

AVG.TFMP 

KFLVIN 

MAX TEMP 
KELVIN 

RMAX 

CMS 

T BULK ( IN) 
KELVIN 

TBULK ( OUT) 
KELVIN 

TWALLM1 

KELVIN 

TWALLM2 

KELVIN 

5.1 

0. 3772F-03 

392.5 

394.8 

0.91727 

366.8 

366.9 

388*6 

387.3 

10.? 

0 . 3777F-03 

39?. 2 

394.4 

0.91747' 

368.4 

368.6 

388.5 

387.5 

15.3 

0.371 7F- 03 

398.? 

400.6 

6.91754 

370.4 

370.7 

393.9 

392.7 

?0. 4 

0.3670F-03 

403.4 

406.1 

0.91761 

372.5 

373.0 

398.8 

397.6 

?5.5 

0. 36? 0F-03 

408.8 

411.7 

0.91770 

374.9 

375.5 

403.8 

402.6 

30.6 

0.3578F-03 

413.8 

416.7 

0.91778 

377.5 

378.2 

408*6 

407.3 

35.7 

0. 3536F-03 

418.5 

421.6 

0.91791 

380.2 

381.1 

413.1 

411.9 

40.8 

0 • 35G2E-0 3 

422. 8 

425.8 

0.91803 

383.0 

384.1 

417*2 

416.1 

45.9 

0.3469F-03 

4?6. 7 

429.8 

0.91820 

385.9 

387.2 

421*0 

420.0 

51.0 

0. 3443F-03 

430.0 

433.0 

0.91838 

388.9 

390.3 

424.3 

423.4 

56.0 

0.3420E-Q3 

432.8 

435.8 

0.91861 

391.8 

393.4 

427*1 

426.4 

61.1 

0. 3406F- 03 

435.0 

438.0 

0.91887 

394.7 

396.4 

429.4 

428.9 

66.? 

0 • 3394F-0 3 

436.6 

439.4 

0.91918 

397.4 

399.4 

431*1 

430.8 

71.3 

0 • 3384F-03 

437.5 

440.1 

0.91955 

400.1. 

402.2 

432.2 

432.2 

76.4 

0*33 84F- 03 

437.8 

440.3 

0.92000 

402.6 

404.9 

432.6 

432.9 

81.5 

0.3384F-03 

437.3 

439.6 

0.92056 

404.9 

407.3 

432.5 

433.0 

86.6 

0. 3394E-03 

436.3 

438.3 

0.92126 

406.9 

409.5 

431.7 

432.6 

91.7 

0. 3409E-03 

434. 5 

436.3 

0.92221 

408.7 

411.4 

430.3 

431.4 

96.8 

0. 34?5F-03 

432. 1 

433.7 

0.92350 

410.2 

412.9 

428.3 

429.8 

101.9 

0. 3454F-03 

428.7 

429.9 

0.92570 

411.4 

414.2 

425.3 

427. 1 

107.0 

0* 3452E- 03 

4?9 . 1 

430.2 

0.92637 

412.5 

415.3 

425.8 

427.7 





are plotted in figure 6 as functions of axial position or elevation. The average helium 3 
density is given in figure 7 as a function of axial position. A comparison of figures 6 
and 7 shows that the average gas density is inversely proportional to the average gas tem- 
perature. At a gas pressure of 41 newtons per square centimeter, the total mass of 
helium 3 in an annulus is 21. 9 milligrams. 


33 


REFERENCES 


1. Davison, Harry W. ; Heath, Colin A. ; and Lowen, Walter: Control of a Nuclear Re- 

actor With a Gaseous Control System. NASA TN D-4093, 1967. 

2. Kreith, Frank: Principles of Heat Transfer. Second ed. , International Textbook Co. , 

1965. 

3. Brown, George G. : Unit Operations. John Wiley and Sons, Inc., 1950. 

4. Shoupp, W. E. : Westinghouse Electric Corporation Periodic Chart. Westinghouse 

Atomic Power Dept. , 1960. 

5. Evans, Robley D. : The Atomic Nucleus. McGraw-Hill Book Co. , Inc. , 1955. 

6. Simmons, Jesse T. : The Physical and Thermodynamic Properties of Helium. Rep. 

No. TR D-9027, Wm. R. Whittaker Co. , Ltd. , July 1, 1957. 

7. Barkas, Walter H. ; and Berger, Martin J. : Tables of Energy Losses and Ranges of 

Heavy Charged Particles. NASA SP-3013, 1964. 


34 



1 



Figure 1. - Program flow sheet for STADDIG. 










































Mass flow rate, 
G i 



Figure 2. - Heat-transfer model. 
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(c) Annular spacing. 

Figure 3. - Geometric configurations. 





Figure 4. - Reactor conditions for sample problem. 
(Drawing not to scale.) 
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Figure 5. - Axial distribution of neutron flux. 
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Figure 6. - Axial temperature distributions in control element 
containing helium 3. 
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Figure 7. - Axial distribution of helium 3 density in annulus. 
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